home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / !RLaB / toolbox / lininsrt < prev    next >
Encoding:
Text File  |  1995-02-19  |  2.1 KB  |  99 lines
  1. % LININSRT Linear insertion between neighboring points.
  2. %    [XO,MASK] = LININSRT(X,N) Inserts N points
  3. %    between each two points of the input vector X.
  4. %    N can be a number or a vector with the length
  5. %    equal to length(X)-1.
  6. %    Returns output vector XO with inserted points 
  7. %    and vector MASK of the same size as XO with
  8. %    0 for old elements and for new ones.
  9.  
  10. %  Kirill Pankratov, kirill@plume.mit.edu
  11. %  01/16/95
  12.  
  13. lininsrt = function ( x , n )
  14. {
  15.   local ( x , n )
  16.  
  17.   % Handle input ...................................
  18.   if (nargs < 2) { n = 1; }
  19.  
  20.   szx = size(x);
  21.   n = n[:];
  22.   ln = length(n);
  23.  
  24.   if (ln == 1) % If scalar, make a uniform vector
  25.   { 
  26.     c = szx[1] == 1;
  27.     c = (!c) * szx[1] + c * szx[2] - 1;
  28.     n = n[ones(c,1)];
  29.     ln = c;
  30.   }
  31.  
  32.   % Check that length of n be equal to the number of
  33.   % columns or number of rows of x minus 1.
  34.   c = (ln == szx - 1);
  35.   if (!any(c))
  36.   {
  37.     error("Second argument length must be size(x,1)-1");
  38.   }
  39.   c = c[2] - c[1];
  40.  
  41.   if (!any(n))  % If no points to be inserted
  42.   {
  43.     xo = x;
  44.     if (c > 0) 
  45.     { 
  46.       mask = zeros(1, szx[2]);
  47.     else
  48.       mask = zeros(szx[1], 1);
  49.     }
  50.     return << xo=xo; mask=mask >>;
  51.   }
  52.  
  53.   % Check if transposition is needed
  54.   is_transpose = 0;
  55.   if (c > 0)
  56.   {
  57.     x = x';
  58.     is_transpose = 1;
  59.   }
  60.   
  61.   szx = size(x);
  62.  
  63.   % Check that all values of n are finite and
  64.   % non-negative.
  65.   mask = find( !finite(n) || n < 0 );
  66.   n[mask] = zeros(size(mask));
  67.  
  68.   % Calculate old and new intervals dx and dxo ........
  69.   dx = diff(x);
  70.   n = n + 1;
  71.   dxo = dx ./ n[;ones(1,szx[2])];
  72.   n = [1; cumsum(n)+1];
  73.   lxo = n[szx[1]];  % Length of the output vector
  74.  
  75.   % Calculate mask .............................
  76.   mask = zeros(lxo, 1);
  77.   mask[n] = ones(size(n));
  78.   a = cumsum(mask);
  79.  
  80.   % Calculate output vector (or matrix) ........
  81.   n = ones(lxo-1,1);
  82.   xo = [x[1;]; x[n;]+cumsum(dxo[a[1:lxo-1];])];
  83.  
  84.   % Make sure "old" values of x are exact
  85.   % (possible round-off error)
  86.   xo[find (mask)] = x;
  87.  
  88.   mask = !mask; % Invert mask (1 - for new points)
  89.  
  90.   % Transpose back if necessary ................
  91.   if (is_transpose)
  92.   {
  93.     xo = xo';
  94.     mask = mask';
  95.   }
  96.  
  97.   return << xo=xo; mask=mask >>;
  98. };
  99.